Smoothing models

Get data - smooth disturbed and control plots together - AGB growth

Code

library(mgcv); library(data.table); library(gratia); library(marginaleffects); library(ggplot2); library(dplyr); library(tidyr); library(zoo); library(slider); library(here)

#####test Paracou####

obs <- read.csv("G:/My Drive/cesab bioforest/Data/aggregated_data_v9.csv")

setDT(obs)

obs=obs[obs$Site=="Paracou"]


setnames(obs, "Plot", "subplot")

obs$plot=gsub("_.", "",obs$subplot )

obs <- obs[!Year%in%1984:1990]

obs <- dcast(
  obs,
  subplot + plot + Year ~ variable,
  value.var = "value"
)


##### upload data 
load("../clim_by_site.rda")

clim_all=rbindlist(lapply(clim_by_site[["Paracou"]], function(x) x$climate[[1]]),
                   idcol = "plot")

#"Paracou", "Mbaiki", "Corinto", "SUAS", "Lesong", "Ulu Muda", "Sungai Lalang", "Tene", "Jari"

baseline <- clim_all[year >= 1981 & year <= 2010]


clim_stats <- baseline[, .(
  mean_tmax = mean(tmax),
  sd_tmax   = sd(tmax),
  mean_vpd   = mean(vpd),
  sd_vpd     = sd(vpd),
  mean_srad  = mean(srad),
  sd_srad    = sd(srad),
  mean_def  = mean(def),
  sd_def    = sd(def)
), by = month]

clim_all <- merge(clim_all, clim_stats, by = "month", all.x = TRUE)

clim_all[, `:=`(
  z_tmax = (tmax - mean_tmax) / sd_tmax,
  z_vpd  = (vpd  - mean_vpd)  / sd_vpd,
  z_srad = (srad - mean_srad) / sd_srad,
  z_def  = (def  - mean_def)  / sd_def
)]

#climate anomalies follow inventories. For Paracou, mean per year before 1996, two years after

clim_all[, year_bin :=
           ifelse(year <= 1995,
                  year,                     # keep original year
                  1997 + ((year - 1996) %/% 2) * 2   # 2-year bins
           )]


census_anom <- clim_all[,.(
  mean_z_tmax = mean(z_tmax, na.rm = TRUE),
  mean_z_vpd = mean(z_vpd, na.rm = TRUE),
  mean_z_srad = mean(z_srad, na.rm = TRUE),
  mean_z_def = mean(z_def, na.rm = TRUE)
),
by = .(plot, year_bin)]

#AGB values before 1991 are not consistent
census_anom=census_anom[year_bin>=1991] 


merged <- merge(
  obs,
  census_anom,
  by.x = c("Year","plot"), by.y=c("year_bin","plot")
)


control=c(1,6,11,13,14,15)

merged$Treatment <- ifelse(merged$plot %in% control, "control", "disturbed")

merged$Treatment <- factor(merged$Treatment)
merged$plot <- factor(merged$plot)
merged$subplot <- factor(merged$subplot)

GAM k=5, k=8, flexible k and bs = “fs”

NAs for agb_growth - plots 13,14 and 15 miss data (?)

Code

merged <- merged %>%
  group_by(subplot) %>%
  mutate(
    fitted_GAM5 = predict(
      gam(agb_growth ~ s(Year, k = 5), data = cur_data()),
      newdata = cur_data()
    )
  ) %>%
  ungroup()


merged <- merged %>%
  group_by(subplot) %>%
  mutate(
    fitted_GAM8 = predict(
      gam(agb_growth ~ s(Year, k = 8), data = cur_data()),
      newdata = cur_data()
    )
  )%>%
  ungroup()

merged <- merged %>%
  group_by(subplot) %>%
  mutate(
    fitted_GAM_flexibleK = predict(
      gam(agb_growth ~ s(Year, k = -1), data = cur_data()), #allows model to choose the best k for each subplot
      newdata = cur_data()
    )
  )%>%
  ungroup()

p1=ggplot(merged, aes(x = Year, y = fitted_GAM5, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="GAM-k=5") + theme(legend.position = "none")


plot(merged$agb_growth, merged$fitted_GAM5)

Code

p1.1 <- ggplot(merged, aes(Year, agb_growth - fitted_GAM5)) +
  geom_point() +
  ggtitle("GAM residuals")

p1 + p1.1

Code

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_GAM5)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("GAM5 residuals vs Tmax")

GAM k=8

Code

p11=ggplot(merged, aes(x = Year, y = fitted_GAM8, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~ subplot, scales = "free_y")+
  
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="GAM-k=8") + theme(legend.position = "none")

p11

Code

plot(merged$agb_growth, merged$fitted_GAM8)

Code

p11.1 <- ggplot(merged, aes(Year, agb_growth - fitted_GAM8)) +
  geom_point() +
  ggtitle("GAM residuals")

 
p11.1

Code

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_GAM8)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("GAM8 residuals vs Tmax")

GAM flexible k

Code

p111=ggplot(merged, aes(x = Year, y = fitted_GAM_flexibleK, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~ subplot, scales = "free_y")+
  
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="GAM-k flexible") + theme(legend.position = "none")

p111

Code

plot(merged$agb_growth, merged$fitted_GAM_flexibleK)

Code

p111.1 <- ggplot(merged, aes(Year, agb_growth - fitted_GAM_flexibleK)) +
  geom_point() +
  ggtitle("GAM residuals")

p111.1

Code

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_GAM_flexibleK)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("GAM8 residuals vs Tmax")

GAM bs = “fs”

Code

#test GAM that is done on the full models (bs = "fs")

merged <- merged %>%
  mutate(
    fitted_GAM8_fs = predict(
      gam(agb_growth ~ s(Year, subplot, bs = "fs", k = 8), data = cur_data()),
      newdata = cur_data()
    )
  )

Code

p111=ggplot(merged, aes(x = Year, y = fitted_GAM8_fs, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  theme_bw() +
  facet_wrap(~ subplot, scales = "free_y")+
  labs(x = "Year", 
       y = "agb_growth", title="GAM-k=8 fs") + theme(legend.position = "none")

p111

Code

plot(merged$agb_growth, merged$fitted_GAM8_fs)

Code

p111.1 <- ggplot(merged, aes(Year, agb_growth - fitted_GAM8_fs)) +
  geom_point() +
  ggtitle("GAM residuals")

p111.1

Code

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_GAM8_fs)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("GAM5 residuals vs Tmax")

Code

#GAM captures early and late patterns faithfully, best edges behaviour. more variance on early years not linked to the model

slider

complete=TRUE (reject incomplete windows, NA at the beggining and the end), window per observation, not true “moving time windows”, window is just row-based.

Code

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide_complete = slide_dbl(agb_growth, mean, .before = 2, .after = 2, .complete=TRUE))%>%
  ungroup()

Code

p3=ggplot(merged, aes(x = Year, y = fitted_slide_complete, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales = "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="slide - .before = 2, .after = 2, window per observation")+ theme(legend.position = "none")

p3

Code

plot(merged$agb_growth, merged$fitted_slide_complete)

Code

p3.1 <- ggplot(merged, aes(Year, agb_growth - fitted_slide_complete)) +
  geom_point() +
  ggtitle("slide residuals")

p3.1

Code

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_slide_complete)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide_complete residuals vs Tmax")

**complete=TRUE, window per Year (NAs can be in the middle - anywhere Year gaps break the window)

Code

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide_completeYear = slide_index_dbl(agb_growth, Year, mean, .before = 2, .after = 2, .complete=TRUE))%>%
  ungroup()

Code

p33=ggplot(merged, aes(x = Year, y = fitted_slide_completeYear, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales= "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="slide - .before = 2, .after = 2, window per year")+ theme(legend.position = "none")

p33

Code

plot(merged$agb_growth, merged$fitted_slide_completeYear)

Code

p33.1 <- ggplot(merged, aes(Year, agb_growth - fitted_slide_completeYear)) +
  geom_point() +
  ggtitle("slide residuals")

p33.1

Code

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_slide_completeYear)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide_completeYear residuals vs Tmax")

complete=FALSE, no NAs, use less points

Code

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide = slide_dbl(agb_growth, mean, .before = 2, .after = 2, .complete=FALSE))%>%
  ungroup()

Code

p333=ggplot(merged, aes(x = Year, y = fitted_slide, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales = "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="slide - .before = 2, .after = 2 - partial windows at edges")+ theme(legend.position = "none")

p333

Code

plot(merged$agb_growth, merged$fitted_slide)

Code

p333.1 <- ggplot(merged, aes(Year, agb_growth - fitted_slide)) +
  geom_point() +
  ggtitle("slide residuals")


p333.1

Code

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_slide)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide residuals vs Tmax")

polynome

use ALL years to fit the curve at each year

Code

merged <- merged %>%
  group_by(subplot) %>%
  mutate(fitted_poly2 = fitted(lm(agb_growth ~ poly(Year, 2), data = cur_data(),na.action = na.exclude)))%>%
  ungroup()

Code

p4=ggplot(merged, aes(x = Year, y = fitted_poly2, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales="free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="Quad poly")+ theme(legend.position = "none")

p4

Code

plot(merged$agb_growth, merged$fitted_poly2)

Code

p4.1 <- ggplot(merged, aes(Year, agb_growth - fitted_poly2)) +
  geom_point() +
  ggtitle("Poly residuals")

p4.1

Code

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_poly2)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("Poly residuals vs Tmax")

root mean square error

Code

rmse_results <- merged %>%
  group_by(subplot) %>%
  summarise(
    RMSE_GAM8   = sqrt(mean((agb_growth - fitted_GAM8)^2, na.rm = TRUE)),
    RMSE_GAM8_fs   = sqrt(mean((agb_growth - fitted_GAM8_fs)^2, na.rm = TRUE)),
    RMSE_GAM_flexiblek   = sqrt(mean((agb_growth - fitted_GAM_flexibleK)^2, na.rm = TRUE)),
    RMSE_GAM5   = sqrt(mean((agb_growth - fitted_GAM5)^2, na.rm = TRUE)),
    RMSE_SLIDE  = sqrt(mean((agb_growth - fitted_slide)^2, na.rm = TRUE)),
    RMSE_SLIDE_COMPLETE  = sqrt(mean((agb_growth - fitted_slide_complete)^2, na.rm = TRUE)),
    RMSE_SLIDE_COMPLETE_YEAR  = sqrt(mean((agb_growth - fitted_slide_completeYear)^2, na.rm = TRUE)),
    RMSE_PolyQuad  = sqrt(mean((agb_growth - fitted_poly2)^2, na.rm = TRUE))
  ) %>% 
  ungroup()

rmse_results <- rmse_results %>%
  mutate(
    best_method = colnames(.)[
      apply(select(., -subplot), 1, function(x) which.min(replace(x, is.na(x), Inf))) + 1
    ]
  )

as.data.table(rmse_results)

    subplot RMSE_GAM8 RMSE_GAM8_fs RMSE_GAM_flexiblek RMSE_GAM5 RMSE_SLIDE
     <fctr>     <num>        <num>              <num>     <num>      <num>
 1:     1_1 0.6602780    0.6579413          0.6604191 0.6597762  0.6862318
 2:     1_2 0.5440815    0.5416930          0.5440815 0.5440815  0.6040639
 3:     1_3 0.6609940    0.6557081          0.6609940 0.6609940  0.7424553
 4:     1_4 0.5800805    0.5803164          0.5800805 0.5800805  0.5960658
 5:    10_1 0.7849051    0.8456235          0.7851315 0.7825946  0.7973029
 6:    10_2 0.6821700    0.6956048          0.6822746 0.6818336  0.6831134
 7:    10_3 0.7672751    0.7549499          0.7673827 0.7682380  0.7282196
 8:    10_4 0.7838802    0.8222049          0.7843187 0.7819798  0.7669645
 9:    11_1 0.7357129    0.7247332          0.7357488 0.7351039  0.7745976
10:    11_2 0.7142709    0.6952320          0.7142709 0.7142709  0.6943462
11:    11_3 0.8031754    0.7794432          0.8031754 0.8031754  0.8220464
12:    11_4 0.8146776    0.7937353          0.8146776 0.8146776  0.7326208
13:    12_1 0.8704537    0.8754618          0.8706103 0.8701286  0.7861002
14:    12_2 0.7417773    0.7371563          0.7417773 0.7417773  0.7541769
15:    12_3 0.9156824    0.9188433          0.9158479 0.9150931  0.8331223
16:    12_4 0.8799091    0.8639448          0.8799091 0.8799091  0.8558761
17:    13_1 0.9668969    0.9638537          0.9671027 0.9692105  0.8783367
18:    13_2 0.9318155    0.9050310          0.9318155 0.9318155  0.8934179
19:    13_3 0.6107574    0.7332620          0.6129402 0.6181148  0.6323810
20:    13_4 0.7650064    0.7648394          0.7648404 0.7651549  0.6167025
21:    14_1 0.7517864    0.7661681          0.7518816 0.7533733  0.7121477
22:    14_2 0.4812247    0.5133870          0.4807792 0.4831050  0.4409424
23:    14_3 0.6529933    0.6871827          0.6524651 0.6535501  0.5654100
24:    14_4 0.6397579    0.6414662          0.6396102 0.6406264  0.5366713
25:    15_1 0.7796617    0.7620071          0.7796617 0.7796617  0.5999784
26:    15_2 0.6667609    0.6523123          0.6668362 0.6667603  0.6222884
27:    15_3 0.6817246    0.6731469          0.6816596 0.6823092  0.7056651
28:    15_4 1.0042997    0.9972733          1.0045377 1.0056790  0.7683992
29:     2_1 0.3835665    0.6019324          0.3921007 0.6227459  0.5444585
30:     2_2 0.4542383    0.5151127          0.4547022 0.4572000  0.4443997
31:     2_3 0.6446780    0.7004875          0.6436109 0.6427212  0.6508805
32:     2_4 0.6203826    0.6101165          0.6203826 0.6203826  0.5947166
33:     3_1 0.7336957    0.7332287          0.7336957 0.7336957  0.6931810
34:     3_2 0.4190214    0.4154055          0.4190214 0.4190214  0.4327896
35:     3_3 0.5668004    0.9629471          0.5964462 0.7456780  0.6978298
36:     3_4 0.4789723    0.4931243          0.4791108 0.4785594  0.4325711
37:     4_1 0.3263275    0.3896101          0.3273419 0.3281933  0.3346153
38:     4_2 0.3738629    0.3762309          0.3737674 0.3736529  0.4240364
39:     4_3 0.6534133    0.7128752          0.6529899 0.6552418  0.6276590
40:     4_4 0.5099975    0.5735599          0.5102925 0.5119338  0.4938520
41:     5_1 0.6145998    0.5986149          0.6145998 0.6145998  0.5774668
42:     5_2 0.4819202    0.4776803          0.4819951 0.4826530  0.4603537
43:     5_3 0.6389234    0.6381846          0.6389234 0.6389234  0.7019902
44:     5_4 0.7110561    0.7032076          0.7110561 0.7110561  0.6971347
45:     6_1 1.0523562    1.0247354          1.0524894 1.0525408  1.0084021
46:     6_2 0.7472778    0.7234551          0.6246955 0.7472778  0.7026307
47:     6_3 0.8233057    0.8158325          0.8233057 0.8233057  0.8410777
48:     6_4 0.7529320    0.7841702          0.7560583 0.7293806  0.7669631
49:     7_1 0.6001234    0.5973120          0.6001234 0.6001234  0.6059278
50:     7_2 0.7373637    0.7359204          0.7373637 0.7373637  0.7936859
51:     7_3 0.6102879    0.6101203          0.6102879 0.6102879  0.6375042
52:     7_4 0.6972089    0.6979671          0.6972089 0.6972089  0.6674518
53:     8_1 0.6704648    0.6765673          0.6701081 0.6632378  0.5724463
54:     8_2 0.7272951    0.7864109          0.7278266 0.7216743  0.7005143
55:     8_3 0.6863311    0.6734833          0.6864467 0.6855547  0.6281695
56:     8_4 0.6313004    0.6696157          0.6306451 0.6302679  0.6508777
57:     9_1 0.6122220    0.6115483          0.6122220 0.6122220  0.5848557
58:     9_2 0.9248327    0.9117891          0.9248327 0.9248327  0.9307054
59:     9_3 0.5965190    0.5913531          0.5964939 0.5966956  0.5728574
60:     9_4 0.7001091    0.6947434          0.7001091 0.7001091  0.6837235
    subplot RMSE_GAM8 RMSE_GAM8_fs RMSE_GAM_flexiblek RMSE_GAM5 RMSE_SLIDE
    RMSE_SLIDE_COMPLETE RMSE_SLIDE_COMPLETE_YEAR RMSE_PolyQuad
                  <num>                    <num>         <num>
 1:           0.6924368                0.6620793     0.6516454
 2:           0.6558556                0.5574081     0.5352876
 3:           0.7302887                0.6020966     0.6561903
 4:           0.6077942                0.5521897     0.5793231
 5:           0.7615287                0.8883787     0.7631782
 6:           0.5688687                0.6474647     0.6665390
 7:           0.6330230                0.7418817     0.7522101
 8:           0.7733975                0.8172911     0.8085605
 9:           0.6619965                0.7246492     0.7102420
10:           0.6526901                0.6654209     0.6958772
11:           0.6402996                0.6816309     0.7745037
12:           0.5623786                0.6018615     0.7902398
13:           0.6936347                0.9333572     0.8510292
14:           0.6736246                0.6983954     0.7400901
15:           0.7329491                1.0766141     0.8982756
16:           0.7499191                0.8705708     0.8643039
17:           0.8533502                0.9836345     0.9410003
18:           0.8904807                0.9827753     0.9016847
19:           0.4058699                0.6067523     0.6737868
20:           0.5866013                0.6967027     0.7444854
21:           0.6788017                0.7558908     0.7276886
22:           0.3884776                0.3321044     0.4667363
23:           0.5026993                0.5559724     0.6330926
24:           0.4515229                0.5204747     0.6219419
25:           0.5157961                0.5710494     0.7528862
26:           0.6637365                0.5446195     0.6430975
27:           0.6263003                0.6668891     0.6579428
28:           0.7643632                0.8526297     0.9706199
29:           0.4984113                0.4239263     0.6122367
30:           0.4632715                0.4124832     0.5254208
31:           0.5823183                0.5416240     0.6288648
32:           0.6498110                0.6627187     0.6044826
33:           0.7493849                0.8052941     0.7236745
34:           0.4018967                0.3380455     0.4172664
35:           0.6317521                0.6565301     0.9045353
36:           0.4409795                0.5441225     0.4836408
37:           0.3718257                0.3777425     0.3868553
38:           0.4424177                0.3714502     0.3632610
39:           0.5018227                0.7070766     0.6354674
40:           0.5429312                0.5322199     0.5234642
41:           0.6168627                0.6321946     0.5946510
42:           0.4977874                0.4988723     0.4721498
43:           0.6431314                0.6411788     0.6383874
44:           0.7800434                0.8189429     0.7015761
45:           0.8832837                0.8802776     1.0070110
46:           0.7221726                0.6071806     0.7162201
47:           0.5719515                0.7388309     0.8130082
48:           0.6417194                0.6385596     0.7735994
49:           0.6295481                0.6167540     0.5991195
50:           0.8030928                0.7717808     0.7367826
51:           0.6110155                0.6849377     0.6102132
52:           0.6693527                0.8425252     0.6969681
53:           0.5828610                0.7311958     0.6814875
54:           0.7502980                0.8807528     0.7636499
55:           0.6672173                0.7487241     0.6749501
56:           0.7189026                0.6477941     0.6291198
57:           0.6333529                0.5900400     0.6117652
58:           0.9564598                0.9876889     0.9120191
59:           0.5418454                0.5428450     0.5772581
60:           0.5723792                0.6445077     0.6963614
    RMSE_SLIDE_COMPLETE RMSE_SLIDE_COMPLETE_YEAR RMSE_PolyQuad
                 best_method
                      <char>
 1:            RMSE_PolyQuad
 2:            RMSE_PolyQuad
 3: RMSE_SLIDE_COMPLETE_YEAR
 4: RMSE_SLIDE_COMPLETE_YEAR
 5:      RMSE_SLIDE_COMPLETE
 6:      RMSE_SLIDE_COMPLETE
 7:      RMSE_SLIDE_COMPLETE
 8:               RMSE_SLIDE
 9:      RMSE_SLIDE_COMPLETE
10:      RMSE_SLIDE_COMPLETE
11:      RMSE_SLIDE_COMPLETE
12:      RMSE_SLIDE_COMPLETE
13:      RMSE_SLIDE_COMPLETE
14:      RMSE_SLIDE_COMPLETE
15:      RMSE_SLIDE_COMPLETE
16:      RMSE_SLIDE_COMPLETE
17:      RMSE_SLIDE_COMPLETE
18:      RMSE_SLIDE_COMPLETE
19:      RMSE_SLIDE_COMPLETE
20:      RMSE_SLIDE_COMPLETE
21:      RMSE_SLIDE_COMPLETE
22: RMSE_SLIDE_COMPLETE_YEAR
23:      RMSE_SLIDE_COMPLETE
24:      RMSE_SLIDE_COMPLETE
25:      RMSE_SLIDE_COMPLETE
26: RMSE_SLIDE_COMPLETE_YEAR
27:      RMSE_SLIDE_COMPLETE
28:      RMSE_SLIDE_COMPLETE
29:                RMSE_GAM8
30: RMSE_SLIDE_COMPLETE_YEAR
31: RMSE_SLIDE_COMPLETE_YEAR
32:               RMSE_SLIDE
33:               RMSE_SLIDE
34: RMSE_SLIDE_COMPLETE_YEAR
35:                RMSE_GAM8
36:               RMSE_SLIDE
37:                RMSE_GAM8
38:            RMSE_PolyQuad
39:      RMSE_SLIDE_COMPLETE
40:               RMSE_SLIDE
41:               RMSE_SLIDE
42:               RMSE_SLIDE
43:             RMSE_GAM8_fs
44:               RMSE_SLIDE
45: RMSE_SLIDE_COMPLETE_YEAR
46: RMSE_SLIDE_COMPLETE_YEAR
47:      RMSE_SLIDE_COMPLETE
48: RMSE_SLIDE_COMPLETE_YEAR
49:             RMSE_GAM8_fs
50:             RMSE_GAM8_fs
51:             RMSE_GAM8_fs
52:               RMSE_SLIDE
53:               RMSE_SLIDE
54:               RMSE_SLIDE
55:               RMSE_SLIDE
56:            RMSE_PolyQuad
57:               RMSE_SLIDE
58:             RMSE_GAM8_fs
59:      RMSE_SLIDE_COMPLETE
60:      RMSE_SLIDE_COMPLETE
                 best_method

lm

Code

vars <- grep("^fitted_", names(merged), value = TRUE)

results <- data.frame(
  method = vars,
  intercept = NA,
  slope = NA,
  p_value = NA
)



for (i in seq_along(vars)) {
  
  v <- vars[i]
  model <- lm(as.formula(paste0(v, " - agb_growth ~ Year")), data = merged)
  s <- summary(model)
  
  results$intercept[i] <- coef(s)[1, 1]
  results$slope[i]     <- coef(s)[2, 1]
  results$p_value[i]   <- coef(s)[2, 4]
}

results

                     method     intercept         slope   p_value
1               fitted_GAM5  2.926468e-11 -1.459840e-14 1.0000000
2               fitted_GAM8  3.746029e-11 -1.375740e-14 1.0000000
3      fitted_GAM_flexibleK -5.365398e-11  2.672661e-14 1.0000000
4            fitted_GAM8_fs -7.298793e-01  3.640826e-04 0.8738790
5     fitted_slide_complete  8.845040e+00 -4.412354e-03 0.1303205
6 fitted_slide_completeYear -5.540808e+00  2.780665e-03 0.3211441
7              fitted_slide  4.140980e+00 -2.064655e-03 0.3527055
8              fitted_poly2 -4.420647e-12  2.205123e-15 1.0000000

temporal signal on residuals

Code

methods <- c(
  "fitted_GAM5",
  "fitted_GAM8",
  "fitted_GAM8_fs",
  "fitted_GAM_flexibleK",
  "fitted_slide",
  "fitted_slide_completeYear",
  "fitted_slide_complete",
  "fitted_poly2"
)

for (m in methods) {
  merged[[paste0("resid_", m)]] <- merged$agb_growth - merged[[m]]
}

pairwise_corr <- function(df, col) {
  wide <- tidyr::pivot_wider(df[, c("subplot","Year", col)],
                             names_from = subplot, values_from = all_of(col))
  M <- as.matrix(wide[, -1])          # remove Year
  C <- cor(M, use = "pairwise.complete.obs")
  mean(C[upper.tri(C)], na.rm = TRUE)
}

temporal_signal <- data.frame(
  method = methods,
  mean_pairwise_corr = sapply(paste0("resid_", methods),
                              function(x) pairwise_corr(merged, x))
)

temporal_signal

                                                   method mean_pairwise_corr
resid_fitted_GAM5                             fitted_GAM5          0.3102782
resid_fitted_GAM8                             fitted_GAM8          0.3110250
resid_fitted_GAM8_fs                       fitted_GAM8_fs          0.3041428
resid_fitted_GAM_flexibleK           fitted_GAM_flexibleK          0.3130019
resid_fitted_slide                           fitted_slide          0.3169307
resid_fitted_slide_completeYear fitted_slide_completeYear          0.2853856
resid_fitted_slide_complete         fitted_slide_complete          0.2960472
resid_fitted_poly2                           fitted_poly2          0.3120015